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C/^ Abstract 

In this paper we propose an approach for simultaneous identification 
of the absorption density and the speed of sound by photoacoustic mea- 
surements. Experimentally our approach can be realized with sliced pho- 
toacoustic experiments. The mathematical model for such an experiment 
is developed and exact reconstruction formulas for both parameters are 
presented. 

g 1 Introduction 

In this paper we propose an approach for simultaneous identification of the 
' ^ absorption density and the speed of sound by photoacoustic measurements. 

*;0 In standard photoacoustic experiments the object is uniformly illuminated by a 

short electromagnetic impulse. Recently sectional photoacoustic imaging tech- 
niques [T71 ini HO] and according reconstructions techniques H] have been 
developed. In sectional experiments thin hyperbola shaped like regions of the 
specimen are imaged (see Figure [T]) sequentially. Experimentally, one can per- 
form sectional photoacoustic imaging by focused illumination combined with 
focusing detectors. The focused illumination is achieved by using cylindrical 
lenses in front of the object. Moreover, contemporary focusing ultrasonic detec- 
tors have a spherical or cylindrical shape, thus the detector surface plays the role 
of the acoustic lens. The generated ultrasonic wave is refracted by a suitable 
acoustic lens such that out-of-plane (center of the hyperbola shaped regions) 
signals are generally weak and can be neglected. Thus essentially only signals 
emerging from the imaging plane are collected at the detector. This justifies 



in 

ON 

o 



X 



*The work has been supported by the Austrian Science Fund (FWF) within the national 
research network Photoacoustic Imaging in Biology and Medicine, project S10505-N20. 

^Karlsruhe Institute of Technology, Department of Mathematics, KaiserstraBe 89, 76128 
Karlsruhe, Germany 

t Computational Science Center, University of Vienna, NordbergstraBe 15, A-1090 Vi- 
enna, Austria, and Johann Radon Institute for Computational and Applied Mathemat- 
ics (RICAM), Austrian Academy of Sciences, Altenbergerstrafie 69, A-4040 Linz, Austria 
(otmar . scherzerOunivie . ac. at) . 



1 



that the illumination can be assumed restricted to a single plane (line in 2D). 
However, such a model requires a low scattering coefficient of the sample (which 
model organism specimens like the Zebra fish have). Opposed to sectional 3D 





Illumination Region 



The object has to be shifted in z-direction 
is therefore illuminated for each horizontal line. 



Fig. 1: Conventional sectional photoacoustic imaging. 



photoacoustic imaging, where stacks of complementary two-dimensional pro- 
jection images are produced, the proposed approach for simultaneous imaging 
consists in performing overlapping sliced imaging by rotation and translation 
of the specimen. This, also generates enough data for reconstructing the two 
independent parameter functions, speed of sound and absorption density, (see 
Figure pi). By now the approach presented here is by far from being experi- 



Fig. 2: Sectional imaging in all directions produces enough data to reconstruct both imaging 
functions. 

mentally economical in the sense that a quick dimension analysis shows that we 
acquire much too many data - in K'^ we reconstruct two 3-dimensional functions 
from eventually six dimensional data. It is the goal of this paper, however, not 
to present the most economical approach, but to show that it is possible to 
derive exact reconstruction formulas for both imaging parameter functions. 

The literature on reconstruction formulas and back-projection algorithms for 
photoacoustic imaging is vast. Wang et al. developed reconstruction formulas 
for cylindrical, spherical, and planar measurement geometries in a series of pa- 
pers [571 [151 [55] ; and recently many more algorithms based on reconstruction 




Illumination Region 




The object has to be shifted and rotated, and 

is therefore illuminated for each hyperbola (hyperboloid). 
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formulas have been developed (see the survey [H]). As it becomes transpar- 
ent below we can make use of reconstruction formulas from photoacoustics- 
essentially we make extensive use of inversion formulas for the spherical mean 
operator. Exact reconstruction formulas for the wavespeed function in ultra- 
sound reflectivity tomography are based on the Born approximation to the wave 
equation, which is valid for moderately varying speed of sound. Reconstruction 
formulas for ultrasound reflectivity tomography have already been derived by 
Norton & Linzer O IS] and are also based on inversion formulas for the spher- 
ical mean operator. The possibility of exact inversion in both fields supports to 
derive exact inversion formulas for both parameters. 

The proposed method is a hybrid and quantitative imaging methods (see [3l [HI 
[5]) for some recent surveys. Most closely related to our approach is the work 
of Stefanov and Uhlmann |24] which presented a photoacoustic experiment for 
recovering either the sound speed or absorption density. A substitute to our 
work is also the recovery of the absorption density function for inhomogeneous 
wave speed p1.fTH[23]. 

The reconstruction formulas for simultaneous imaging utilize techniques from 
reflectivity imaging and photoacoustic imaging. In the current state of research 
we can provide exact reconstruction formulas, but when we try to be economical, 
i.e., using less slicing experiments, we would require reconstruction formulas, 
which have not been developed so far. This is highlighted in Section [2j The 
outline of this paper is as follows: In Section [2] we discuss the mathematical 
model and discuss the non-economicality of proposed model, which leaves room 
for further improvement of the results. In Section [3] we transform the model 
into the Fourier domain, from which in Sections |3.2[ |3.1| exact reconstruction 
formulas are derived in 2D, 3D, respectively. The Appendix provides the exact 
definitions used in this paper. 

2 Model 

In M", we consider the following Cauchy problem for the wave equation: 

4 dttu - Au = in M" x M>o , 

it(x,0) = f{x)5r,B{x) inM", 
dtu{x,Q) = inM". 

where c = c{x) denotes the speed of sound and 5rfi{x) = 5{d\st{x,E{r,9))) 
where E{r^ 9) is the [n — l)-dimensional hyperplane with distance r from the 
origin and orientation 9 e 5"^^. This is the photoacoustic equation with sliced 
illumination in the plane E{r,9). It is the ultimate aim to reconstruct c and / 
from measurements of u on some hypersurface T (see below). 

We consider the Born approximation; that is, we expand u formally with respect 
to the contrast function q:= Ijc? — 1 and consider only terms of order at most 
q. This leads to the decomposition {t w u + w where u = u^'^ is the solution of 
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the wave equation 



dttu - Aw = in M" x M>o , 

u(a;,0) = f{x)5r,e{x) in E" , 

9tu(a;,0) = in M" , 

and V — v^'^ solves 

dttv — Av = —q{x) dttu in M" x M>o , 
v{x,Q) = inE", 
dtv{x,0) = inE". 

The dependence of v on r,9 is through the source function u. 
It is our goal to reconstruct q and / from measurements of 

m'--''{x,t) = u'-''^{x,t) + v''-\x,t), e r X (0,r), (2.1) 

where F is a (n — l)-dimensional hypersurface in E". Currently everything is 
fixed to complete measurements on the whole surface T. 

We make the following assumptions: 

1. r = dB for some open and connected set B C E" with smooth boundary. 
As typically in photoacoustics, we consider B to be a sphere, a circle, or 
a half plane. 

2. / is the sum of some known initial distribution /o with compact support 
and some unknown term fi to be determined. 

3. The supports of fi and q are both contained in Q for some bounded domain 

ncB, 

4- /o, /i and q are smooth and f{x) = fo{x) + fi{x) ^ for all x E B. 



Dimensionality Analysis 



In the ?i-dimensional setting, we record measurements (cf. (2.1)) for every 
{x,t) G r X (0,T) and every (r,6') G (0,cx.) x 5*"-^ That is, the recorded 
data are 2n-dimensional. The data to be recovered are /i and q, which are two 
n-dimensional functions. In general we think that we record too many data 
for the purpose of reconstructing /i and q - however, since we rely on Radon 
transforms techniques for exact inversion, we accept this disadvantage for the 
mathematical studies. 



3 Fourier Reconstruction Formulas 

In this section we derive exact reconstruction formulas for /i and g. First we 
derive a general inversion formula, which is then evaluated differently in 2D and 
3D. 
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Fig. 3: Schematic representation of the domains and supports. 



In the following we omit the superscripts (r, 9) for the sake of convenience of 
notation. Let 



u{x, k) = J^u{x, k) 



u{x,t)e"'Ut, a;GM", A; G . 



be the Fourier transform - note that u vanishes for t < 0. We assume that u is 
bounded in M" x M>o such that u{x, •) G L^(0, oo) for every x G M". Then we 
note that u has the following properties (by the theorem of Palcy- Wiener) : 

(a) u{x, •) has a holomorphic extension into C+ := {z & C : > 0} for all 

xgW, 

(b) for every a; G M" it holds that J^^\u{x,ki +ik2) — u{x,ki)\^dki — > as 
k2 tends to zero, 

(c) u{-, k) is bounded in M" for all k G C+. 
We have 



u{x, k) = Tu{x, k) 



u{x,t) e'^^dt 



1 

ifc 



u{x,t) e^'^^lZ. 



-lu{x,0) - ^ 



dtu{x,t) e'^'dt 
dtu{x,t) e'''*^°° - 



U=o 



POO 

/ dttu{x,t)e'^*dt 
Jo 



= "^"(2;,0) - ^Au{x,k), 



or in other words 



Au[x,k) + k^u{x,k) = iku{x,0). 



(3.1) 
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Moreover, we have 



v{x,k) = Fv{x,k) = --p- / vtt{x,t) e"''* dt 



[Av{x, t) — q{x) dttu{x, t)] dt 



_ 1 

= -^M{x,k) + ^ J^°° dtMx,t)e"''dt 
= -^Av(x,k) + ^Au{x,k), 

and thus 

Av{x,k) + k'^v{x,k) = q{x)Au{x,k) 

~ —k'^q{x)u{x,k) + ik q{x) u{x , 0) . 

Let be the (radiating) fundamental solution of the Helmholtz equation Aii + 
k'^u = in M" for k & C with 5z > 0; that is, in particular for n = 2,3, 

^Hl^'^\k\x-y\) forn = 2, 

where H^^^ denotes the Hankel function of the first kind and order zero. Then 
one particular solution of (3.1 1 is given by 

u{x,k) = -ifc / u{y,0)^k{x,y)dy 



= -ifc / fiy)'^k{x,y)ds{y). 

Jy(£E{r,d) 

We note that this particular solution satisfies the properties (a), (b), (c) from 
the beginning of this section. Any other solution is of the form u{-, k) + w{-, k) 
where w satisfies Aw + k^w — in M". The requirement that the solution 
satisfies the properties (a), (b), (c) from the beginning of this section implies 
that w vanishes. Indeed, for fc S C+ the function w{-, k) has to be bounded by 
property (c), therefore its Fourier transform (in the distributional sense) with 
respect to x e M" satisfies 

F.^w{y, k) [fc2 - = for all y e M" 

which implies Fxw{-,k) ~ because fc^ — does not vanish for fc e C+. 
Therefore, also w{-, fc) = in M" for all fc e C+ and, by the continuity property 
(b), also w(-,k) = in M" for all fc e C with 3? > 0. 
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We also have 

v{x, k) 

= I q{y)u{y,t)^k{x,y)dy - ifc / q{y) f{y) 5r,e{y)'^k{x,y) dy 
= / q{y) / f{z)^kiy,z)ds{z)<^>k{x,y)dy 

■/R" J z£E(rfl) 

-ik q{z) f{z)<^k{.y,z)ds{z) 

Jz(EE(r,0) 



= -ifc Jf{z) k^ J q{y)^k{y,z)^k{x,y)dy + g(z)$fe(x,z) ds{z) . 

zeE{r,e) 

In summary, we have 

m'^'^{x,k) = vJ''^[x,k) + v'^-'^{x,k) = u{x,k) + v{x,k) 

= -ifc/ f{z)\k^ q{y)^k{y,z)^k{x,y)dy 

JzeE(r,e) L JR" 

+ {q{z) + l)<^k{x,z) ds{z) 

= R[{f{-)L{x,;k))\{r,e), 

where R[f]{r,9) is the (n — l)-dimensional Radon transform of / in direction 
(r, 9) and 

L{x,z,k) = i-'^kf qiy)<l>k{y,z)<^>kiy,x)dy - ik{q{z) + 1) <i>k{x, z) (3.2) 

for .T e r, z e s, fc e M, x 7^ z. 

In other words, we have 



J'^'^{R^\u + v)){x,z,t) = f{z)L[x,z,t) , 



(3.3) 



for a; e r, z e S, fc e M, x 7^ z. 

Therefore, from the knowledge of rh^'^{x, k) — u{x, k) + v{x, k) for all x € F, 
fc G M, r > 0, and 6 e S*"^^ we can determine 

/(z) L(x, z, t) = (/o(z) + /i(z)) L{x, z, 

for all a; G r, z € B with x ^ z, and t > 0. 

3.1 3D Domain 

For n = 3 we recall that 

"^kix^y) 



gifc|x-j(| 

4:T:\x — y\ 



, x^y- 



Thus from (3.2 1 it follows that 
1 



L{x,z,k) = (-ifc)^ 



167r2 



^ik{\y-z\ + \x-y\) 

\x-y\\y- z\ 



■dy - ifc(g(z) + l) 



^ik\x-z\ 

At:\x — z\ 



Taking the inverse Fourier transform with respect to k gives 
327r3 \Jr3 \x - y\\y - z\ J _^ 



8772 |x- 

Twice Integration with respect to t gives 

/ / L{x, z,t) drdsi 
Jo Jo 

- -H{s2-\x-z\) 



87r2|a; — z 

I \ 

'All) 

2>1'k'^\z — y\\y — X 



f 

/ H\y-z\ + \y-x\-S2) ^^ 31 ""''ii ^ds{y) 

JM3 ... 



where we used the Heavyside function H{t) = 1 for t > and H{t) = for 
r < 0. The initial conditions for the second term vanish since for t = the 
domain of integration (w.r.t. y) is empty. 

Integrating once more, wo get 

L{x,z,T)dTdsids2 (3.4) 



f H{s2-\x-z\)ds2+J\f[q]ix,z,t). 
X- z\ Jo 



q{z) + 1 
87r2| 



Now, we note that 

/ H{s2-\x- z\)ds2 = {t-\x-z\)+ := {t-\x-z\)H{t-\x-z\), 
Jo 

and thus 

nx,z,t)= ^['^^\ {t-\x-z\)+ +Af[q]{x,z,t). (3.5) 

on \x — z\ 

Both terms vanish for t < |a; — ^;|. 

We recall that M{x, z, t) := f{z)^{x, z, t) is known from measurement data for 

all a; G r, z G B, and t > 0. First we take z e B\Q. Then q{z) = and 
f{z) = fo{z), thus '^{x,z,t) — M{x,z,t)/ fo{z) is known and 

fo{z) 8w^\x-z\ 
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Again, the left hand side is known. Now we let z tend to x. The limit on the 
right hand side exists, thus also on the left hand side, and therefore for all a; e F 
and t > we have 



lim 



M{x,z,t) 
fo{z) 



t 



8t:^\x~z\ 



1 



d{\y-x\-t/2)qiy)dsiy) 
M2[q]{x,t/2) , (3.6) 



1 



where A^2['7] is the spherical mean operator (see Appendix). 
Thus the reconstruction algorithm is as follows: 



1. Calculate the product f{z)L(x,z,t) from (3.3) for all xGT, 
z G B and t>0. By integration (see (3.4)), this yields the 
knowledge of M{x,z,t) :— f (z)'^ {x , z , t) for all x G F, z Cz B and 
t > 0. 



2. Solve (3.6) for g in by inverting the spherical mean opera- 
tor . 



3. Compute ^(a;,z,i) for all a; G F, z £ fl and t>0 from (3.5). 

4. Finally, compute / from f{z) = M(x,z,t)/'^(x,z,t) for zGfi. 

Remark 3.1 The operator Af[q]{x, z,t) is the rotational ellipsoidal mean oper- 
ator with focal points x and z. Thus the integral equation {3.5); that is, 



■^{x,z,t) 



qjz) + 1 
8n^\x~ z\ 



{t-\x-z\)+ + U[q]{x,z,t)., 



can, for instance, be considered as a fixed point equation for q involving the 
ellipsoidal mean operator and can be solved by the fixed point iteration 

^{x,z,t) = |4r^(i- |a:-z|)+ + Af[qn-i]{x,z,t),n=l,2,... 

Ellipsoidal mean operators have been studied in John's book US^ . 



3.2 2D Domain 

In we recall that the fundamental solution is given by 



'^k{x,y) 



H^o'\k\x - y\) , x^y. 



We have to compute the inverse Fourier transform of the product $fe(x, y) ^kiz, y) 
(see (3.2|) and use the convolution theorem. First we have that 



1 



T-'{H^o\-\^-y\))it) = <! 2n^t^-\x-y\^ 

0, 



, t> \x-y\^ 
t<\x-y\ 
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The convolution for the inverse Fourier transform is given by the formula 
^'^Ua) = f *9; that is, 



^-'ifgm - = / f{T)g{t~T)dT. 

J — oo 

Therefore, 

(^^) ^-\H^,'\-\x-y\)H^,'\-\z-y\))it) 

i2^I\x-l\ '^^2_|^_j^|2 7(t_r)2-|^-j;|2 t-\y - z\> \x-y\, 

0, i - - z| < |x- - yl , 



and thus from (3.2 1 

L{x,z,t) (3.7) 

for < — |y — z| > |x — yl where H denotes again the Heaviside function. Again, 
as in 3D, both terms vanish for t < \x ~ z\. 

Letting z tend to x (thus t > \x — z\ for z sufficiently close to x) yields 
L{x, X, t) 

I t-\v-z\ \ 

\2\y-x\<t \x-y\ / 



Note that q{x) vanishes for a; G F. Using polar coordinates and denoting 

: dr , t > 2r > , 



1 1 



k{t,r) := <^ Jr V^^^ ^{t-rY- 



2 



0, 0<i<2r, 



we therefore get 



*/2 /.2. / /cospNN /-'/^ 



k{t,r)r I q\x + r\ . dpdr = 27r / A^i [g] (x, r) fc(t, r) r dr 

Jo V Vsmp// 7o 

with the one dimensional circular mean operator A^i. Using the notations 
k{t,r) :— rk{2t,r) and 

L{x,x,T)dTdsids2 -2Tr{t\og{t)-t) , (3.8) 
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we get, for fixed x, the Volterra integral equation for Mi[q]{x, ■): 



= / k{t,r)Mi[q]{x,r)dr, te{0,T]. 
Jo 



(3.9) 



We express fc as a complete elliptic integral: We introduce the variable (j) S 
I + (I — J') sincj). Then dr = — cos (j)d<j) and, with 



(-7r/2,7r/2] by r 
t' = t/2, 



2 2 

r — r 



(t — r)(r + r) 
= [(t' - r) + {t' ~ r) sin </>] [(t' + r) + {t' - r) sin 
{t-Tf-r^ = {t-T -r){t-T + r) 

= [{t' - r) - {t' - r) sin </«] [(t' + r) - (f' - r) sin ( 

For the product we conclude that 

[r^ - r^] [{t rf - r^] 
= [{f - rf - {t' - rf sin^ 0] [{t' + rf - {t' - rf sin^ cf\ 
= [t' - rf cos" cj) [{t' + rf - [t' - rf sin^ 0] . 



Therefore, 



k{t,r) 



7r/2 



^/^ ^{f + rf - [f - rf sin^ (j) 
d<j) 



r/2 



and thus 



and where 



k{t,r) — pK{l — p) with p 



K{a) 



■it/2 



\/l — sir? (j) 



2r 
t + r 



\a\ < 1, 



denotes the complete elliptic integral. 

Let 6 > such that 6 < dist(r, fi). Since Aii[q]{x,T) = for r < 5 we can 



consider the integral equation (3.9) on [S,T]; that is 



^{x,t) = I k{t,r)Mi[q]ix,r)dr , t e [S,T] 
s 



(3.10) 



Now we note that k{t,t) = K{0) = 7r/2 and k e C^iA) where A = {(i,r) G 
[S,Tf : r < t}. Therefore, wc can transform the equation of the first kind into 
one of the second kind by differentiating (3.10). This yields 



$(x,t) = ^Mi[q]{x,t) + I ^^k{t,r)Mi[qKx,r)dr, te[S,T]. (3.11) 
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Volterra integral equations of the second kind with smooth kernel and non- 
vanishing diagonal have a unique solution, see, e.g., [131 Section 3.3]. 

After calculating the spherical mean operator i [g] , it can be inverted by us- 
ing standard inversion formulas, which exist for various geometries: For B — 
B]j{0) C M^, analytical reconstruction formulas have been derived by Finch, 
Haltmeier, Rakesh [6j. For a general domain Q, Kunyansky reduced in |16j the 
reconstruction problem to the determination of the eigenvalues and normal- 
ized eigenfunctions Uk, ||ufe||2 = 1, of the Dirichlet Laplacian — A on f2 with 
zero boundary conditions. Recently Palamadov [U] presented a general ap- 
proach leading to reconstruction algorithms for various geometries. For more 
details see also the Appendix. 

Thus the reconstruction algorithm is as follows: 



1. Calculate the product M(x, z,t) := f{z)L{x, z,t) from (3.3) for all 



xeF, z e B, t > 0. Compute ^{x.t) from (3.8) for all x e T 
and t > 0. 

2. Solve ( |3. for calculating Mi[q]. 

3. Invert the spherical mean operator. 

4 . From ( |3.7[ ) determine L by differentiation. 

5. Finally use f(z)~M{x,z,t)/L{x,z,t) to reconstruct /(z) for zG 

n. 
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4 Appendix 

Notation 

B^{x) denotes the ball of radius R and center x in M". 
Bessel Functions 

The Bcsscl functions, Neumann functions and Hankel functions of order zero 
are defined as 

1 r 



^ 2. 



[ e'" dr forzeC, and 

J — TT 

w 



Yoiz) = ^[^ + ln(z/2)] Jo(z) - (^)'' for z e C with 3 > , 



H^^'^ = Jo + iYo; 
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where 7 denotes Euler's constant, ai = —7, and a„ = —J + '^^^i ^ ^ for n > 2 



From this we conclude that Hq^\~x) = Hq^\x) for a; G M. 

The Fourier Transform 

In this paper we use the foUowing notation: 

/oo 
fit) e"=*dt, 
-00 

denotes the Fourier transform of / and 

fit) := T-\fit) := — / /(fc)e-"'=dfc, 

the inverse Fourier transform. For / G L^(M) or / from the Schwarz space 
S they are defined in the classical sense, for / £ L^(M) they are defined by 
extension using the Parseval formula, and for tempered distributions they are 
defined by duality. 

Examples: 

/•C30 

dik) = ISik) = / e"=*(5(t)dt = e"=° = 1, 



00 



1 f°° 1 
Sit) = J'-^Sit) = — / e-'^*5ik)dk = — 
27r . _^ 27r 



This shows that 



1 1 

lik) = 27r5(fc) and !(<) = — / 6""=* d/fc = 

271" ./_oo 27r 



Furthermore, we have (see pU] ) 

00 

li7W(fca) = f^^0!Ldt forfc>0. 



From iJp"'^^ = iJp^'' (x) for all a; G M we conclude that this formula holds for 
all fc G M with k ^ 0. The right hand side is a Fourier transform. Therefore, 

1 



V'* / [ 0, t<a. 

Radon Transform 

In M", Eir,6) denotes the (n — 1)— dimensional hyperplane with orientation 
6 G and distance r from origin. 

Let / : M" — > i? with support in fl, then the (n— l)-dimensional Radon transform 
of / is defined by 

Rr.eif) = R[f]ir..O) = / f{x)dsix) . 

JE(r,e) 

Thus R[f] is a function from M x S"-'^ into M. 
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Spherical Mean Operator 

In M" the spherical mean operator is defined as follows (see e.g. [5]): 

Mn^i[u]{x,r) / u{x + ry) ds{y) ior x e W\r > 



The spherical mean operator can be written as 

1 



Mn-i[u]{x,r) ^ — _ _ / S{\x - z\ - r)u{z) dz for X £K" ,r > 



r 

In particular 

M2[u]{x,r) = I u{x + ry) ds{y) = \ [ S{\x - z\ - r)u{z) dz , 

A^i[u](a;, r) = — / u{x + ry) ds{tj) = / 5{\x — z\ — r)u{z) dz . 

2tt J si ' 27rr Jr2 

There exist a variety of reconstruction formulas for averages over circular and 
spherical means: 

in 2D 



For averaged data on the sphere of Radius i?, that is on dD. = dBj^{0) C 
E^, analytical reconstruction formulas have been derived by Finch, Halt- 
meier, Rakesh ^ and read as follows 

/(C) = ^A^ J^'' riM,[f]Mr)\og\r^ - \^ - e\^\drds{e)^ 

and 

1 



fiO = 7r / idrrdrMi[f]){0,r)log\r'-\^-0\'\drds{0). 
^i" Jsi Jo 

• For a general domain 17, Kunyansky reduced in [IB] the reconstruction 
problem to the determination of the eigenvalues and normalized eigen- 
functions Ufc, ||ufc||2 = 1, of the Dirichlet Laplacian —A on with zero 
boundary conditions: 

Aufc(0 + AfcUfc(0 = 0, ^en, 

Indeed, if (C, 77) 1— ?> {\£,~v\) is a free-space rotationally invariant Green's 
function of the Helmholtz equation and n{^) denotes the outer unit normal 
vector of dH. at ^ S dil, then 

00 

fc=0 

where 

Mk^ f I rMi[f]{ri,r)Gx,{r){Vuk{v),n{r]))drds{ri). 
Jan Jo 
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in 3D 



We refer to the survey of Finch and Rakesh [8] (see also [7]). In this article 
three reconstruction formulas are documented, which are 

If du {t\M,[f])M)\ 
2t:Ro JoBj^ \x-xo\ 

If" (5ti9tt(A^2[/]))(xo,t)|,^|._,„| 

/W = -7^—^ / 1 1 M^o) , 

^T^Ro JdBl |a;-a;o| 
/(a;) = --^A(/ \x-xo\M2[f]){xo,\x-xo\)ds{xa)] . 



2ttRo 

Here x € B^^^^ and / is compactly supported in this set 
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